NETLOGO_PATH <- '/home/rick2/Apps/netlogo-5.0.5'
MODEL_PATH <- 'viralload.nlogo'
NSim <- 10
NTicks <- 650  # FIXME: tem que ser repetido dentro de simfun()

tic <- Sys.time()

library(parallel)
processors <- detectCores()
cl <- makeCluster(processors)


# initialization
presim <- function(dummy, gui, nl.path, model.path) {
print('presim')
	library(RNetLogo)
	NLStart(nl.path, nl.version=5, gui=gui)
	NLLoadModel(model.path)
}

# the simulation
simfun <- function(x) {
	NTicks <- 650  # FIXME: pq?!
	NLCommand('setup')
	d <- data.frame(T=numeric(NTicks), A=numeric(NTicks), D=numeric(NTicks), A0=numeric(NTicks))
	for(t in 1:NTicks) {
		T <- NLReport('count patches with [pcolor = T]')
		A <- NLReport('count patches with [pcolor = A1 or pcolor = A2]')
		D <- NLReport('count patches with [pcolor = D]')
		A0 <- NLReport('count patches with [pcolor = A0]')
		d[t,] <- c(T,A,D,A0)
		NLCommand('update')
	}
	d
}

# quit
postsim <- function(x)
	NLQuit()


nl.path <- Sys.getenv('NETLOGO_PATH', NETLOGO_PATH)
model.path <- paste(getwd(), MODEL_PATH, sep='/')
parLapply(cl, 1:processors, presim, gui=FALSE, nl.path=nl.path, model.path=model.path)

result <- t(parSapply(cl, 1:NSim, simfun))

parLapply(cl, 1:processors, postsim)
stopCluster(cl)

toc <- Sys.time()
print(toc-tic)

## Process results

M <- apply(result, 2, function(r) apply(matrix(unlist(r), nrow=NTicks), 1, mean))
V <- apply(result, 2, function(r) apply(matrix(unlist(r), nrow=NTicks), 1, sd))
#MIN <- apply(result, 2, function(r) apply(matrix(unlist(r), nrow=NTicks), 1, min))
#MAX <- apply(result, 2, function(r) apply(matrix(unlist(r), nrow=NTicks), 1, max))

# standard error = standard deviation / sqrt(NSim)


library(reshape2)
X <- melt(data.frame(t=1:NTicks, M), id.vars='t')
Y1 <- melt(data.frame(t=1:NTicks, V), id.vars='t')
Z <- data.frame(X, ymin=(X$value-Y1$value), ymax=(X$value+Y1$value))

# change time scales so that the first weeks are displayed with more promeninence
Z$t[Z$t > 12] <- Z$t[Z$t > 12] + 100
Z$t[Z$t <= 12] <- (Z$t[Z$t <= 12]-1) * 100/11

output <- function(filename, width)
	pdf(paste(filename, '.pdf', sep=''),
		paper='special', width=8, height=8/width)

library(ggplot2)
library(gridExtra)

output('virus', 2)
p <- ggplot(Z, aes(x=t, y=value, group=variable, ymin=ymin, ymax=ymax)) +
	geom_ribbon(aes(fill=variable), alpha=0.2) +
	geom_line(aes(color=variable)) +
	xlab("") + ylab(expression(paste('# cells x ', 10^5, sep=''))) +
	scale_color_manual("Cell Type", values=c('blue','green','red','orange')) +
	scale_fill_manual("Cell Type", values=c('blue','green','red','orange')) +
	scale_x_continuous(breaks=c(1,50,100,seq(100+365/7,750,by=365/7)), labels=c(1,6,12,seq(2,13,by=1))) +
	scale_y_continuous(breaks=seq(0,5e5,by=1e5), labels=seq(0,5,by=1)) +
	geom_vline(xintercept=c(100), linetype="dotted") +
	geom_text(data=data.frame(label=c("weeks","years"), x=c(50,400)), aes(label=label, x=x, y=y, group=NULL, ymin=NULL, ymax=NULL), y=-Inf, vjust=3.2, size=4) +
	ggtitle("Cell count dynamics - no treatment")
gt <- ggplot_gtable(ggplot_build(p))
gt$layout$clip[gt$layout$name=="panel"] <- "off"
grid.draw(gt)
dev.off()

